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^: 

a^ ■ ABSTRACT 

rH ' It is weU known that the standard transport equations violate causality when gradi- 

i_^ ' ents are large or when temporal variations are rapid. We derive a modified set of transport 

Tj- I equations that satisfy causality. These equations are obtained from the underlying Boltz- 

,__! I mann equation. We use a simple model for particle collisions which enables us to derive 

*^ ', moment equations non-perturbatively, i.e. without making the usual assumption that the 

O ■ distribution function deviates only slightly from its equilibrium value. We also retain time 

,-H ■ derivatives of various moments and choose closure relations so that the final set of equa- 

S ' tions are causal. We apply the model to two problems: particle diffusion and viscous 

^ ■ transport. In both cases we show that signals propagate at a finite speed and therefore 

"^ ' that the formalism obeys causality. When spatial gradients or temporal variations are 

' I small, our theory for particle diffusion and viscous flows reduces to the usual diffusion and 

^ ! Navier-Stokes equations respectively. However, in the opposite limit of strong gradients 

c^ . the theory produces causal results with finite transport fluxes, whereas the standard theory 

> ■ gives results that are physically unacceptable. We find that when the velocity gradient is 

^ ■ large on the scale of a mean free path, the viscous shear stress is suppressed relative to 

c^ ■ the prediction of the standard diffusion approximation. The shear stress reaches a maxi- 

mum at a finite value of the shear amplitude and then decreases as the velocity gradient 
increases. The decrease of the stress in the limit of large shear is qualitatively different 
from the case of scalar particle diffusion where the diffusive flux asymptotes to a constant 
value in the limit of large density gradient. In the case of a steady Keplerian accretion 
disk with hydrodynamic turbulent viscosity, the stress-limit translates to an upper bound 
on the Shakura-Sunyaev oj-parameter, namely a < 0.07. The limit on a is much stronger 
in narrow boundary layers where the velocity shear is larger than in a Keplerian disk. 
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1. INTRODUCTION 

According to the standard theory of diffusion the particle flux increases hnearly with 
the gradient of the particle density. Thus, the theory predicts that a sharp gradient will 
result in a divergent flux. This violates causality, since the particle flux is obviously limited 
by the flnite particle speed. The standard theory runs into a similar problem when the 
particle density changes on a very short time scale locally (e.g. Morse & Feshbach 1953, 
Narayan 1992). 

A similar limitation is well-known in viscous interactions. In this case the viscous 
stress tensor is linearly proportional to the velocity gradient in the diffusion approxima- 
tion. When this relation is used in the Navier-Stokes equation one finds an instantaneous 
propagation of viscous signals, in violation of causality. Since the particle speed does not 
enter explicitly in the underlying theory, the standard equations cannot be modified in a 
straightforward way to limit the signal propagation speed. 

The diffusion equation and the Navier-Stokes equation are valid only when particles 
have suffered many collisions and their distribution has relaxed to have weak spatial gra- 
dients and slow temporal variations. However, there are physical situations both in the 
laboratory and in astrophysical gases where gradients are large on the scale of a coUisional 
mean free path or temporal changes are rapid relative to the mean collision time. Exam- 
ples include radiation hydrodynamics in optically thin media (Levermore and Pomraning 
1981), electron heat transport in laser produced plasma (Max 1981), and viscous angular 
momentum transport in boundary layers of accretion disks (Popham & Narayan 1992). 
One needs to go beyond the standard fiuid equations to model such conditions. 

Specific prescriptions have been proposed to incorporate causality in individual prob- 
lems. In radiation hydrodynamics, Levermore and Pomraning (1981) introduced a flux- 
limiter such that, regardless of the gradient, the radiative flux is never permitted to exceed 
the product of the radiation energy density and the speed of light. In the accretion disk 
problem, Narayan (1992) proposed a modifying factor for the viscosity coefficient, such 
that in steady state fiows the viscosity vanishes when the fiow speed exceeds the maximum 
random speed of the particles. 

While the above prescriptions have worked well in their individual applications, it 
would be useful to derive a general formalism that automatically yields a causal limit to dif- 
ferent transport phenomena. We present such a formalism in this paper. We base our work 
on the Boltzmann equation which is strictly causal. We make a simple non-perturbative 
approximation to the scattering term in the Boltzmann equation, take successive moments, 
and use appropriate closure relations. The standard diffusion approximation is then re- 
covered if we neglect certain terms involving time derivatives. This approximation is valid 



if temporal variations are slow compared to the collision time of particles, and spatial 
gradients are small, but it breaks down whenever there are rapid variations and the flux 
is limited by causality. When we retain all terms in the moment equations, we obtain a 
causal set of equations. 

The technical discussion in the paper is divided into two main sections. In §2 we 
discuss the diffusion of particles in a fixed background, while in §3 we consider the stress 
tensor. 

The diffusion problem provides a simple test-bed for our approximations since it in- 
volves both temporal and spatial variations of the gas properties. We introduce the basic 
features of our approach in §2.1 where we discuss particle diffusion in one dimension. There 
are two main simplifications which permit us to obtain a set of closed causal equations 
in this case. First, we write the scattering term of the Boltzmann equation in a simple 
non-perturbative form. This explicit form allows us to take moments of the Boltzmann 
equation and to obtain equations that are valid for arbitrarily strong gradients in the par- 
ticle density. This particular approximation of the scattering term is used throughout the 
paper and is a key ingredient of our approach. Second, we close the equations by assuming 
that the second velocity moment of the distribution function is constant. In §2.2 we explore 
the properties of our one-dimensional diffusion equation and verify the validity of these 
assumptions by comparing results with numerical tests. We then generalize the theory to 
three-dimensional diffusion in §§2.3 and 2.4. 

Section 3 describes the viscous transport of momentum in the presence of strong 
velocity gradients and presents the main results of this paper. In §3.1 we derive a causal 
equation for the evolution of the stress tensor. We use the same approximation to the 
scattering term that was employed in §2. However, since the elements of the stress tensor 
are themselves second moments of the distribution function, we close the equations at 
the level of the third moments rather than the second moments. We use the simplest 
approximation allowed, assuming that all third moments vanish. The causal equations we 
thus derive reduce to the Navier-Stokes equation whenever the velocity gradients are weak, 
but can also be used when the gradients are large. We apply the new equations to a number 
of problems involving a steady state shear (§3.2-§3.5), and show that in the presence of 
a large velocity gradient the viscous stress is suppressed compared to the prediction of 
the standard diffusion approximation. We consider the effect of steady advection on the 
viscous shear stress in §3.4 and discuss bulk viscosity in §3.5. In §3.6 we extend our analysis 
to rotating flows. The presence of a Coriolis force introduces the epicyclic frequency into 
the problem, and this leads to a generalization of our formula for the modified shear stress. 
In §3.7 we discuss the implications of this formula for accretion disks and show that the a. 



parameter introduced by Shakura and Sunyaev (1973) has a strict upper limit whenever 
the viscosity is mediated by hydrodynamic interactions. FinaUy, we summarize the main 
conclusions in §4. 

Parts of this paper overlap previous work in the subject. The causal particle diffusion 
equations which we derive in equations (2.1.11) and (2.3.6) have appeared several times 
in the literature (Israel and Stewart 1980 and references therein, Schweizer 1984). Many 
of the previous discussions have been somewhat phenomenological whereas we derive the 
equations through a systematic procedure which reveals clearly the specific assumptions 
which we make. In particular, we discuss the limits of the diffusion theory and identify 
exactly which phenomena are described well by the theory and which are not. Our discus- 
sion of the shear problem, especially the effects of strong shear and advection, appears to 
be largely new. Some of the results on rotating flows have been derived independently by 
Kato and Inagaki (1993) whose preprint we received at a late stage of our work. 



2. PARTICLE DIFFUSION 

2.1 The One-Dimensional Problem 

We introduce our notation and explain our basic ideas and approximations by dis- 
cussing first a one-dimensional problem. Consider a gas of light particles diffusing in a 
fixed background of much heavier particles. Let the light particles be described by a dis- 
tribution function /(t, x, v) where t is time, x is particle position, and v is particle velocity. 
We define various moments of / in the usual way: 

n = fdv, V = — vfdv, v'^ = — v'^fdv. (2.1.1) 

J n J n J 

We assume that the particles experience no accelerations in between scatterings. 
The distribution function / satisfies a Boltzmann equation of the form 

| + _|(„;) = r..„., + /„ (2.1.2) 

where Tscatt describes the effect of scattering and fait, x, v) is the rate at which particles 
are introduced or removed from a phase space element. Let us write Vgcatt as 

Tscatt = r+-r- = -(/o-/). (2.L3) 

T 

The term F" = —f/r represents the rate at which particles are removed from a phase 
space element, where t{v) is the mean free time which is in general a function of velocity. 
r+ = /o/r describes the phase space distribution of the particles after they are scattered. 
Provided we make /o a function of / and normalize /o so as to conserve particles, equation 
(2.1.3) will always be true. Note that this equation does not imply that / is close to 
/o in any way. Indeed many of the cases we consider in this paper involve highly non- 
linear situations where / is very different from /q. In this respect we differ from the usual 
approach to transport theory where one expands / around an equilibrium /o and assumes 
that the deviations are small (see e.g. Krook's equation [Shu 1992]). 

In general, the scattering is a complicated function of velocity which makes the Boltz- 
mann equation very difficult to handle. A major simplification is achieved if we make two 
approximations (see Appendix B of Grossman, Narayan, and Arnett 1993, also Kato and 
Inagaki 1993). First, we assume that r is a constant, independent of velocity. Second, we 
make some simplifying assumptions regarding the post-scattering distribution function /o, 
viz. we assume that each scattering is elastic in the frame of the fixed background, that 



it conserves particles, and that it leads to a total randomization of the initial velocities. 
This allows us to write simple relations for the moments of /q: 

fodv = / fdv = n, / vfodv = 0, / v'^fodv = nv'^. (2-1-4) 

In the same spirit, we assume that the source function fg has zero mean velocity, and write 

fsdv = s{t,x), - vfsdv = 0, - v'^fsdv = v'^. (2.1.5) 



With the above assumptions, we now take the first two moments of equation (2.1.2). 
Integrating equation (2.1.2) over v, we obtain 

- + -(„.)^. (2.1.6) 

Multiplying equation (2.1.2) by v and integrating over v, we obtain 

d _ d 1 _ 

-—(nv) + 7— (nt;2) = nv. (2-1.7) 

Ot ox T 

Equations (2.1.6) and (2.1.7) are two coupled equations describing the evolution of the 
particle moments. Unfortunately, these equations involve three different moments, n, v 
and v'^. Therefore, we need a closure relation. We make the simplest assumption possible, 
namely that the velocity dispersion of the particles, t;^, is a constant, independent of x 
and t: 

v^ = a"^ = constant. (2.1.8) 

This condition is exactly valid if all particles, including those created by the source function 
/s, have a single speed, v = ±cr, and provided the scattering is elastic. In a more general 
situation, this condition may not be satisfied but one might hope that it will perhaps still 
capture the essential features of particle transport phenomena. We discuss later the degree 
to which the approximation (2.1.8) is valid and where it breaks down. Equation (2.1.7) 
now becomes 

«(""' + " & = """• (2-i-»> 

Equations (2.1.6) and (2.1.9) provide a closed pair of equations for the two moments n and 

V. 

If we ignore the time derivative in equation (2.1.9) we get the usual formulation of 
diffusion, where the diffusive fiux nv is given by —Ddn/dx, with a diffusion constant 
D = a'^T. In this approximation n satisfies the standard diffusion equation 

|-.0^. 0..^r. (.00) 



As is well-known, this equation violates causality (e.g. Morse & Feshbach 1953, Narayan 
1992). 

If, on the other hand, we retain the time derivative in (2.1.9) we obtain a modified 
diffusion equation which does preserve causality. Differentiating equation (2.1.9) with 
respect to x and combining with equation (2.1.6) we obtain 



dn 2 f d^n 1 d'^n 



ds 



(2.1.11) 



This equation has been derived previously in the literature using the theory of "transient 
thermodynamics" (Israel and Stewart 1980, Schweizer 1984). The most interesting feature 
of equation (2.1.11) is the presence of the wave operator on the left hand side, which 
ensures that signals cannot propagate faster than the r.m.s. particle speed a. 

The Green's function Gi{t,x) of the one-dimensional diffusion equation (2.1.11) is 
obtained by solving the equation 



dGi 2 d'Gi 

(J T 



dt 



1 d^Gi 



m^r±m 



5{x). 



(2.1.12) 



Morse and Feshbach (1953) have given the solution when the right hand side consists only 
of the first term in the square brackets. Modifying their solution for our case (Schweizer 
1984, Nagel and Meszaros 1985), we obtain 



Gi(t, x) 



-^\h{u) + ^^-^\ + \{5{x-at)+6{x + at)} 



-t/2i 



\x\ < at. 



(2.1.13) 



= 
where Iq and /i are modified Bessel functions and 

„2\ 1/2 



\x\ > at, 



1 



u=— t'-^ 
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(2.1.14) 



Note that the Green's function cuts off for \x\ > at. This demonstrates explicitly that 
equation (2.1.11) satisfies causality with a maximum propagation speed of a. 

Equation (2.1.11) and the Green's function (2.1.13), which are derived under several 
approximations, represent the physical situation exactly in the particular case when all 
particles have the same speed a and the same mean free time r. The 6 functions and the 
causal fronts aX x = ±aT in equation (2.1.13) are of course a consequence of the mono-speed 
assumption. The generalization of the Green's function for an arbitrary distribution of 



speeds, /(Ifl), and an arbitrary dependence of the mean free time, r(|f |), is straightforward 
and requires merely integrating equation (2.1.13) over the particular distributions*. 

One effect that is not described by equations (2.1.11) and (2.1.13) is the smoothing 
of density inhomogeneities by phase mixing (see e.g. Binney & Tremaine 1987). Particles 
with different velocities travel different distances within a collision period and therefore 
wash out inhomogeneities. The damping of inhomogeneities by free streaming of particles 
is not described by equation (2.1.11) since it corresponds to particles with a single speed. 
However, for elastic collisions this problem can be easily fixed by integrating the Green's 
function (2.1.13) over the velocity distribution of particles. 

Another way to improve the description of particle diffusion is to drop the assumption 
of a constant v^ (which crudely speaking corresponds to an isothermal system) in the mo- 
ment equations. We can treat v^ as an independent moment and solve for it by considering 
higher moments of the Boltzmann equation. This extension of the theory describes both 
particle and heat diffusion. We consider up to third moments of the Boltzmann equation, 
and use a closure on the fourth moments of the form v^ = C(^^)^ with ( equal to some 
constant (e.g. C = 3 for a Gaussian distribution). Let us generalize equations (2.1.4) and 
(2.1.5) suitably for the higher moments, 

v^fodv = 0, I v^fdv = Cn(y^f, (2.1.15) 

- / v^fsdv = 0. (2.1.16) 

We then obtain 

— (n^) + — (n^) = s^, (2.1.17) 

ot ox 

^(n^) + C-^[nf^f] = n^. (2.1.18) 

ot ox T 

By taking the x— derivative of equation (2.1.18) and substituting equation (2.1.17) one 
gets a diffusion equation for nv'^, 

'*"^' r(ci;[(n^Pl-SM)=r^ + W). (2.1,19) 



dt V^^ ' ' df^' 'j dt 

which is coupled to the diffusion equation for n 



dn I d'^lnv'^] d'^n\ ds ,„ , „^, 



* The conditions, however, are different in more complex gases, like plasmas. Colli- 
sionless Landau damping may dominate over coUisional dissipation whenever the plasma 
properties undergo strong temporal or spatial variations. 



Note that the hmiting speed for the propagation of a pressure perturbation is greater by 
a factor ~ y^ than that of a density perturbation. 

Equations (2.1.19) and (2.1.20) provide a more accurate representation of diffusion 
than the simpler version of the theory presented earher. Furthermore, these equations 
are causal as is evident from the wave operator on the left side. In fact, these equations 
also describe thermal diffusion since the quantity nv^ in equations (2.1.17) and (2.1.18) 
represents the flux of heat. We do not discuss these extended equations further in this 
paper but concentrate instead on the properties of the simpler equation (2.1.11). 

2.2 Evaluating the Accuracy of the Approximate Theory 

Figure 1 shows the Green's function Gi at two different times, t = O.lr, and lOr. For 
comparison, two other Green's functions are also shown for each case. One is the Green's 
function Gs of the standard diffusion equation (2.1.10). The other is the Green's function 
Gm for a Maxwellian distribution of particles with an r.m.s. one-dimensional velocity a. 
This is calculated by integrating equation (2.1.13) over a Maxwellian in a, and assuming 
r to be independent of v. At the detailed level of the shape of the Green's function, 
neither equation (2.1.10) nor equation (2.1.11) does a particularly good job of fltting the 
Maxwellian Green's function, particularly at early time. This is not surprising since the 
two theories have only one parameter {D) and two parameters (r, a^) respectively, and 
therefore they cannot describe the exact spatial distribution of particles with a continuous 
velocity distribution. 

Although the shape of the Green's function Gi differs signiflcantly from Gm, the 
spatial extent of Gi is close to that of the Maxwellian Green's function Gm at all times 
(see flgure 1). This is in contrast to the Green's function Gs of equation (2.1.10) which 
is much too wide at early times (the usual signature of acausal behavior). We now show 
that the mean square width {x^) of Gi is, in fact, exactly equal to {x^) of Gm- 

Let us deflne {x^) for Gi as follows: 



{x ) = dxx Gi(t, x), where / dxGi{t,x) = 1. (2.2.1) 

Multiplying equation (2.1.12) by a;^ and integrating over dx we flnd 

d{x^) d^{x^) , 

^^ + r^^ = 2ar, (2.2.2) 

whose solution with the appropriate boundary conditions is 

(x^) = 2a^Tt + 2(TV2[exp(-t/r) - 1]. (2.2.3) 

9 



At early times, i.e. for t -^ t, we see that (x^) = a'^t'^ which corresponds to particles 
streaming freely with a speed a. At late times, however, we have (x^) -^ 2a'^Tt which 
corresponds to the usual diffusion limit. 

Since the Green's function Gi corresponds exactly to the case of particles with a 
single speed cr, equation (2.2.3) is the exact result for the evolution of (x^) for a mono- 
speed population of particles. The interesting point is that the expression for (x^) is 
directly proportional to a'^. Therefore, any distribution of velocities which has a mean 
square velocity equal to a^ will evolve exactly according to equation (2.2.3). This equation 
is therefore valid more generally than for just mono-speed particles — all it requires is that 
all the particles should have the same mean free time. 

If particles with different speeds do not have the same mean free time, then obviously 
we do not have a perfect correspondence between (2.2.3) and the exact result for (x^). 
However, even in this case we find that the theory performs quite well. As a demonstration 
of this result, suppose we have particles with a Maxwellian distribution of one-dimensional 
speeds, 

f{v)dv = J^exp(-^]dv, v>0, (2.2.4) 



TTO" 



2a2 



and let us assume that all the particles have the same mean free path /. The mean free 
time is then a function of v, 

t{v)=1/v. (2.2.5) 

To obtain the evolution of {x'^) for this problem, we replace (j"^ by v"^ in equation (2.2.3) 
and integrate over the Maxwellian f{v)dv. This gives 

(^') = 2yf ^/t + 2/2 [exp (^^^ erfc (^^^ - l] . (2.2.6) 

At early and late times this has the following asymptotic dependences: (x^) -^ a'^t^ , t <^ t; 
{x^) — » 2(2/7r)^/2(T/t, t ^ r. We can fit these dependences with the single-speed result 
(2.2.3) provided we choose the mean free time r to be 

r = \[^-l. (2.2.7) 

V TT a 

Figure 2 shows a comparison between the exact result (2.2.6) and the approximate result 
(2.2.3) for this particular choice of r. The agreement is very good. In contrast, the standard 
diffusion equation (2.2.10), which gives {x^) = 2a'^Tt for all t, clearly makes a large error 
for t < T. We expect our theory to give similar good agreement with exact results for 

10 



other distribution functions f{v) or prescriptions for the mean free time t{v), provided we 
define the mean free time r appropriately. 

We thus conclude that the causal diffusion theory developed here provides a good 
representation of particle diffusion, and in particular it predicts the mean square distance 
traveled by particles very well at all times. 

2.3 Particle Diffusion in Three Dimensions 

We work in a Cartesian representation with position indicated by the coordinates 
r = (xi, X2, xs) and velocity hj v = {vi, V2, fa). In the absence of any forces, but with 
scattering, the Boltzmann equation gives 

% + ^(^^^) = ^-«" + A = ^(/o - /) + fs. (2.3.1) 

We use the summation convention on repeated indices. Following the one dimensional 
example of §2.1, we have written the scattering term simply in terms of a mean free time 
r. We also make similar assumptions as in equations (2.1.4) and (2.1.5), viz. 

fodv = n, / vJodv = 0, / fgdv = s, - Vifsdv = 0. 

Taking the zeroth and Vjth moments of (2.3.1) we obtain 

dn d , , 

— + --nv-, = s, 2.3.2 

Ot OXi 

di^^"^^ ^ a^^^^^^ " --nl^J. (2.3.3) 

These are not a closed set of equations because they involve the six second moments ViVj . 
To close the equations, we simplify the second moments by assuming that the velocities 
are isotropically distributed and that the mean square velocity is independent of t and r : 

viVj = -a'^Sij, o"^ = constant. (2.3.4) 

This approximation is technically inconsistent, particularly at early times. While the mean 

square velocity itself can be constant, for instance in the case of a mono-speed population of 

particles, the assumption of isotropy is rather drastic and we may expect some pathological 

behavior in the solutions in certain limits. Nevertheless, the resulting theory is simple and 

causal, and therefore we explore it further. Substituting equation (2.3.4) into equation 

(2.3.3) we find 

d a^ dn 1 

11 



Combining equations (2.3.2) and (2.3.5) we then obtain 

dn 1 9 / 9 3 d'^n\ ds 

_-_,V(v^n--j^j=. + .g^. (2,3.6) 

As before, we find a wave operator appearing in the equation which enforces causahty. 
Rather surprisingly, the relevant wave speed appears to be cr/y^- This is however some- 
what deceptive and the true speed turns out to be a in certain situations as we show in 
§2.4. For reference we note that the standard diffusion equation in three dimensions is 

^-Dy''n = s, D = -a\. (2.3.7) 

Morse and Feshbach (1953) show that the three-dimensional Green's function 6*3 (t, r) 

of the operator in equation (2.3.6) can be derived from Giit^x) discussed in §2.1. The 

relation is 

G3(t,r) = --^-^Gi(t,r), (2.3.8) 

Zixr dr 

with a in equation (2.1.13) replaced by a /\/?> because of the modified wave operator 
in equation (2.3.6). The appearance of derivatives of the 5-function in G^ means that 
this Green's function can give unphysical negative particle densities under certain circum- 
stances. This may be somewhat surprising considering the fact that the one dimensional 
Green's function was not unphysical in any way. There is, however, a significant physical 
difference between the one-dimensional and three-dimensional problems. In one dimen- 
sion, we showed that equation (2.1.11) corresponds to a particular real physical situation, 
namely the case of mono-speed particles. This is unfortunately no longer true in three 
dimensions. There is no physical three dimensional system which behaves exactly as de- 
scribed by equation (2.3.6). In particular, a mono-speed population of particles does not 
obey this equation. The reason is associated with the closure assumption that the sec- 
ond moment tensor of the velocity is isotropic (eq. 2.3.4). This assumption is obviously 
violated close to an impulsive point source before collisions take place. 

As a final comment we note that if we consider a particular case when all the gradients 
are restricted to one direction, say the xi axis, then the three-dimensional equation (2.3.6) 

reduces to 

dn 1 9 fd n 3 d n\ ds . , 

This is identical to the one-dimensional equation (2.1.11) except for the identification 
o^ -^ a"^/?). As we have seen in §2.1 and §2.2, the one-dimensional problem behaves 
quite reasonably (it has no negative densities), and therefore the three-dimensional version 
derived here should be well behaved when gradients are nearly uni-directional. 

12 



In §3, we derive a causal equation of evolution for the shear stress and consider a 
number of flows with gradients in the velocities. By analogy with the particle diffusion 
problem we expect well behaved results in the shear flow problem when gradients are in the 
same direction everywhere. All the situations we consider in §3 do have gradients limited 
to a single direction. 

2.4 Accuracy of the Three-Dimensional Diffusion Model 

We have seen that the Green's function G3 contains the derivative of a 5-function 
which yields an unphysical density at early times. Thus, the solutions to equation (2.3.6) 
cannot be accurate in detail. Following the discussion in §2.2, we can however ask whether 
equation (2.3.6) is accurate in some more limited sense. For instance, does it fit some 
spatial moment of the particle distribution? The answer is that equation (2.3.6) does 
indeed do an excellent job of predicting the evolution of the mean square distance (r^) 
traveled by particles. 

To show this we set s = S{t)S^{f) in equation (2.3.6), multiply the equation by r^, 
and integrate over Airr'^dr. We find that 

+ r^^T^=2(7V, (2.4.1 



dt dt^ 



whose solution for the given boundary conditions is (cf. eq. 2.2.3) 

(r^) = 2a^Tt + 2(TV2[exp(-t/r) - 1]. (2.4.2) 

This result shows that (r^) varies as a'^t'^ at early time (t <^ r), corresponding to free 
streaming at a speed a. (This is a little surprising since equation (2.3.6) has an apparent 
speed limit of a j \f?> in the wave operator. The Green's function Gj,{t,r) too cuts off for 
r > atj \f?). Mathematically, it is the derivative of the 5- function in G3 which leads to an 
r.m.s. particle distance of crt even though the Green's function cuts off at a smaller radius.) 
At late time, (r^) has the standard dependence due to diffusion, viz. (r^) = 2cr^rt. 

Motivated by the early time behavior of (r^) in equation (2.4.2), let us find the en- 
semble averaged evolution of (r^) for a set of particles with a mean square velocity a^ 
which are injected at the origin at a time t = 0. The Brownian motion of each particle is 
described by Langevin's equation (Reif 1965), 

dv V -^, . . , X 

- = — +F{t), (2.4.3) 

at T 

The first term on the right hand side represents a frictional force, where r here is the same 
as the mean free time in our theory and is taken to be a constant as per our assumptions. 

13 



The second term, F{t), is a rapidly fluctuating force with zero mean. We assume that this 
force is arranged so that the mean square velocity (f^) = a"^ is independent of time (again 
as in our model of the scattering). Multiplying equation (2.4.3) by r, averaging over the 
ensemble of particles, and noting that {f ■ F) = 0, we get 

{^(f-v))=a'--{f-v), (2.4.4) 

at T 

With the identity (f ■ v) = |((i(r^)), equation (2.4.4) admits the solution, 

(r2) = 2a^Tt + 2(7^2 [exp(-t/r) - 1]. (2.4.5) 

which is identical to equation (2.4.2). This means that equation (2.3.6) predicts the evo- 
lution of the mean square distance (r^) exactly for any population of particles with mean 
square velocity a^ and constant mean free time r. This result clarifies in what sense equa- 
tion (2.3.6) is a good representation of diffusion in three dimensions — it is not very good 
at predicting the detailed shape of the density distribution but it does provide an exact fit 
of the mean square distance diffused by particles. 

As an extension, we ask what would happen if the mean free times of the particles are 
not all equal. As before we consider the particular case when all particles have the same 
mean free path /, i.e. 

t{v) = l/v, (2.4.6) 

and where the distribution of speeds is Maxwellian in three dimensions: 

47r / 3 \^/^ / 3^2 \ 

f[v)dv = — { — \ v^e^v{--^\dv, v>Q. (2.4.7) 

After replacing a by v, one can integrate equation (2.4.2) over this f{v) to find 




(r^)^,,I^^.Jt-,i^± 



exp(w^)erfc(ty) 
a 



-2/2, w = ^^. {2A1 

>i /v6a 



This result has the following asymptotic behaviors: (r^) -^ a^ t^ for t <^ t; and (r^) — * 
2(8/37r)^/2(TJ^^/t for t ^ t. Comparing with equation (2.4.2) we see that we can obtain 
agreement in the behavior of (r^) in the two limits if we choose the mean free time r in 
equation (2.4.2) to be 

r = ./f-L. (2.4.9) 

Figure 3 shows a comparison between the exact result (2.4.8) and the approximate result 
(2.4.2) when r is defined as in (2.4.9). The agreement is excellent, showing that our 
diffusion theory predicts (r^) very well even when the mean free time is not constant. 
Also shown in figure 3 is the result from the standard diffusion equation, viz. (r^) = 
2(8/37r)^/2a"J^^/t, which is very inaccurate at early time. 
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3. STRESS TENSOR 

3.1 Evolution Equation for the Stress Tensor 

Having introduced our approach through the relatively simple problem of particles 
diffusing through a fixed background, we now consider the case when particles scatter off 
one another. We decompose the velocity of a particle Vi into a mean velocity vl plus a 
relative velocity Ui: 

Vi=Wi + Ui. (3.1.1) 

Particles experience an acceleration a = (ai, a2, 03), where the a^ may be functions of time, 
position, and in general velocity as well. We ignore the possibility of particles being added 
or removed from the system. We thus have the Boltzmann equation 

We depart from the discussion of §2 in the properties we assume for the scattering. 
Because the particles scatter off one another, the relevant frame in which we should specify 
the properties of Tgcatt is not some fixed background frame but rather the local rest frame 
of the gas of particles. Since each scattering event conserves momentum, the mean velocity 
of the particles must be preserved. Thus 

(W)o = - / Vifocfv = v~, (117)0 = - Uifod^v = 0. (3.1.3) 

n J n J 

Furthermore, we specify the second moments of the post-scattering distribution function 
/o to be 

-If 2 

{u~uj)o = - UiUjfod^v = {l-^) — Sij, a'^ = vqWi = ul + ul+ul. (3.1.4) 

In essence we assume that the scattering completely isotropizes the particle velocities so 
that immediately after scattering, (i) there are no off-diagonal velocity correlations such as 
{uiU2)o, and (ii) the mean square velocity is independent of direction, i.e. {u\)o = (^3)0 = 
(wDo- A new feature in equation (3.1.4) is the introduction of the parameter ^, which 
allows for the possibility of inelastic scattering (^ = corresponds to elastic scattering). 
The point is that particles heat up when there are stresses present, and in real gases some 
of this heat is lost from the system by, for example, radiative processes. We introduce ^ to 
model this kind of "cooling" which has to be modeled if we are interested in considering 
steady state situations. 
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Let us take the zeroth and first moments of tlie Boltzmann equation (3.1.2). Tfiese 
give tfie standard continuity and momentum equations: 

| + |-(„W)=0. (3,1.5) 

dv7 dv7 ^ 9 , , , , 

— f- + Vi—-^ = ttj -—[nuiUj). (3.1.6) 

ot oxi n oxi 

In equation (3.1.6), there are two kinds of acceleration, a mean acceleration oj due to 
the external forces and a contribution from the gradient of the stress tensor nuiUj. In 
order to be general and to permit velocity-dependent external forces, we Taylor expand 
the acceleration aj in the form 

«j=a7+-^«^. (3-1-7) 

where aJ is the mean force and the second term is the excess force on a particle due to 
its relative velocity. Equation (3.1.7) is quite general and the only approximation made is 
that the expansion in velocity has been truncated at the linear term. 

We now derive equations for the evolution of the stress tensor. To do this it is conve- 
nient to switch from the distribution function /(t, {xi}, {vi}) to another "relative" distri- 
bution function /i?(t, {xi}, {ui}) (cf. Kato 1970, Grossman et al. 1993) which is a function 
of the relative velocities Ui. This distribution function satisfies a Boltzmann equation 

where iii is the rate of change of the relative velocity of a particle as we follow it along its 

trajectory: 

^ _ dtti dWi dWi 

Ui = Vi - Vi = tti + ^—Uj - ^- - [Vj + Uj)- — 
OVj ot OXj 

= -^—Uj - - — Uj + -——[nuiUj). (3.1.9) 

OVj OXj n OXj 

We now substitute equation (3.1.9) into equation (3.1.8) and take the UjUk moment to 
obtain the equations of evolution of the stress components, 

d d d ( dv ' dci ■ x 

Qrinujuj:) + -Q—inWi uju^) + —{nv~v~uj:) = - ( — ^ - — ^ j nujPn~ 

- I -^ -K— I nUjUi H na djk nUjUk- (3.1.10) 
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Equivalently the equations for second moments are obtained by combining (3.1.5) and 
(3.1.10), 

9 , , _ d . , 1 d , , (din daj\ 

at oxi n oxi V 9^i 9vi J 






As usual, the above equations involve higher order moments (in this case the third 
moments UiUjUk) and we need additional relations to close the set of equations. We 
faced such requirements for closures in the discussion on diffusion in §2. There we needed 
closures for the second moments and we made the simplest assumption possible, viz. that 
the second moments are constant and isotropic. Since the third moments are odd in 
velocity, the simplest assumption here is to set all the third moments to zero: 

UiUjUk = 0. (3.1.12) 

This is a reasonable approximation whenever the gradients of the second moments are 
small. More precisely, we expect odd moments in the velocity to be negligible whenever 
there is sufficient spatial symmetry in the neighborhood of the fluid element under consid- 
eration. Note that equation (3.1.12) does not require velocity gradients to be small. Indeed 
much of what we discuss later in the paper deals with large velocity gradients. We do how- 
ever make sure that there is sufficient "right-left" spatial symmetry in the model problems 
we solve so that equation (3.1.12) will be valid. In the Appendix we derive a conditions on 
the second and third derivatives of the mean velocity which are necessary to make relation 
(3.2.11) viable for a steady shear flow. We should mention that there are other possible 
closures for the third moments. For instance, one could use a "diffusion approximation" 
to relate the third moments to gradients of the second moments, e.g. uf oc ufrdul/dxi, 
and so on. Relations like this are reasonable, but they ultimately correspond to setting 
time derivatives to zero in the third moment evolution equations and using a Gaussian-like 
closure on fourth moments (see equation [2.1.18]). The problem with the neglect of time 
derivatives is that it may lead to a violation of causality. As an analogy, note that in 
the particle diffusion problem, if we neglect the time derivative in (2.1.9) then we obtain 
the usual diffusion approximation which gives the acausal equation (2.1.10). In the same 
way we expect a diffusion-like closure relation on third moments to lead to a violation of 
causality. 

Equations (3.1.5), (3.1.6) and (3.1.11) with third moments set to zero provide 10 equa- 
tions (one continuity equation, three momentum equations, and six stress equations) that 
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describe the evolution of the 10 moments: n, ni, and WjUj. The Navier-Stokes equation 
is a special case of these equations. When the spatial and temporal derivatives are small 
on a collision scale (3.1.11) gives the standard linear relation between the stress and the 
velocity gradient. This result can in turn be substituted into equation (3.1.6) to get the 
Navier-Stokes equation. While the Navier-Stokes equation does not satisfy causality, our 
full equations can handle strong velocity gradients and satisfy causality. The combined 
equations (3.1.6) and (3.1.11) can potentially replace the Navier-Stokes equation in situa- 
tions where the gradients or the temporal variations are arbitrarily large. In the following 
sections we explore the causal properties of these equations especially for the particular 
case of steady flows with a one-dimensional velocity gradient. More work, however, is 
needed in order to examine the full regime of applicability of equations (3.1.5), (3.1.6) and 
(3.1.11) under general conditions. Note that, at the level of our approximation there is 
no independent equation of state. Our system of particles behaves essentially like an ideal 
gas, and whatever prescription we may use for the cooling will determine the degree to 
which the gas is non-adiabatic. Deviations from an ideal gas behavior could be introduced 
by allowing for additional degrees of freedom in the particles but we do not explore this 
possibility here. 

It is useful to write the stress tensor explicitly as the sum of the diagonal terms and 
the off-diagonal terms. The diagonal terms can be decomposed into two parts: 

nuf =p + n[uf --a j , ^ - ^' (3.1.13) 

where the flrst term is the isotropic pressure. The second term describes the degree of 
anisotropy in the diagonal elements and gives rise to bulk viscosity. The off-diagonal 
components of the stress tensor constitute the shear stress. We discuss the physics of the 
shear and bulk viscous stresses within our formalism in the following sections by considering 
various limiting cases of equation (3.1.11). 

Before we proceed to discuss steady state conditions let us show that the stress equa- 
tion (3.1.11) satisfles causality. For simplicity, consider a flow, in the absence of external 
forces, where the mean velocity at each point is oriented along the X2 axis and its magni- 
tude 1)2 is a function of xi. Since this flow is divergence-free we take the particle density n 
to be constant. In addition we take uf ~ cr^/3 to be constant, which is a valid assumption 
for a sufficiently weak shear {uiU2 ^ c'^)- Under these conditions the (1,2) component 
of equation (3.1.11) together with the X2-component of the momentum equation (3.1.6) 

yield, 

dv2 f d'^v^ 3 d'^v^ 



at "» ( ifef - JJIJF I = "• <^-i-"> 
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where z/q = a^r/3 is the viscosity coefficient. Equation (3.1.14) has the same structure 
and shares the same causal properties as the diffusion equation (2.3.11). Thus, based on 
the discussion in §2 we expect our shear flow theory to be causal and well behaved in one 
dimension. 

3.2 Steady State Linear Shear 

We begin our discussion of steady state shear flows by considering a flow where there 
are no external accelerations, i.e. a^ = 0, and where the mean velocity at each point is 
oriented along the X2-axis. We assume that the magnitude of the velocity is independent 
of X2 and xs, and varies linearly with xi, i.e. 

IJ7 = 0, v^ = 2Axi, W = 0. (3.2.1) 

We assume that the system is in a steady state, and take n to be a constant, as is consis- 
tent with the divergence- free nature of the flow. Under these conditions, the momentum 
equation (3.1.6) shows that the various second moments have no spatial gradients. Fur- 
thermore, by the symmetries of the problem, the only non- vanishing second moments are 
uf, w|, w|, UiU2- Thus, the four non- vanishing components of the second moment equation 
(3.1.11) are given by 

(1 - 0(7^ - 3wf = 0, (3.2.2) 

{1 - C)(J^ - ^vi - I^AtWu^ = 0, (3.2.3) 

(1-^(7^-3^ = 0, (3.2.4) 

2ATwf + wTU^ = 0. (3.2.5) 

We add to these four equations the deflnition of a^ given in equation (3.1.4). 
If we sum equations (3.2.2) — (3.2.4) we obtain the following relation, 

^ 1 

-uTu2 ■ 2A = - ■ -a^ . (3.2.6) 

T 2 

The left-hand side of this equation is the work done (stress x shear) per particle per unit 
time by the shear stress. The right hand side represents the rate at which kinetic energy 
is lost per particle as a result of cooling. These two quantities should obviously be equal 
in a steady state and equation (3.2.6) shows that our system of equations is physically 
consistent. More importantly, it reveals the need for the cooling parameter ^, since without 
its presence there are no steady state solutions to the equations (except for the trivial case 
of a zero shear) . 
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Going back to equations (3.2.2)-(3.2.5), we can solve for ^ to obtain 
and from this we solve for all the second moments in terms of a^: 

?=i±^-^ (3.2.9) 

Note that in deriving these results we did not restrict the linear shear parameter 2A in 
any way and therefore these formulae are valid even for a large linear shear as long as the 
spatial derivatives of the shear are sufficiently small (see the Appendix). However, we did 
assume a steady state. Therefore, the results may not be applicable when the flow varies 
on a collision time. Also, we assumed that all particles have the same mean free time. 
Even when this is not the case, we expect that we can choose r appropriately, so that 
our equations will describe the behavior of the system well. This point was demonstrated 
for the particle diffusion problem in §2.2 and §2.4, and will be discussed further in the 
next subsection. Similarly, there is no problem in principle with spatial variations of r. 
Finally, we note that the stress tensor depends on the form of the cooling function assumed. 
However, the main qualitative features of the results, such as the decrease in the stress 
U1U2 at large shear amplitude, appear to be universal. 

Equations (3.2.7)-(3.2.10) show that the nature of the steady state stress tensor de- 
pends on the magnitude of the dimensionless parameter At. In the limit of a weak shear. 
At <^ 1, we find that the velocity distribution of the particles is essentially isotropic, i.e. 
-uf ~ wl ~ w| ~ cr^/3, and the shear stress is given by 

a'^T dv2 

nuiU2 = —n ■ 2A = — nz/Q ■ 7; — • (3.2.11) 

3 oxi 

This is the standard linear relation between the viscous stress and the velocity shear. The 
kinematic coefficient of viscosity is 

z/0 = ^. (3.2.12) 

The Navier-Stokes equation is based on a relation of the form (3.2.11) for the shear stress. 
This relation, plus an equivalent one for the bulk viscosity which we discuss in §3.5, are 
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introduced into the momentum equation to give a closed set of dynamical equations for 
the flow. We see that this approximation is valid only when the shear is weak, i.e. only 
when the shear frequency 2A is small compared to the scattering frequency 1/r. For larger 
velocity gradients we have to use the more complete set of equations presented here. 

When the shear is strong, i.e. At ^ 1, we find a completely different behavior than 
the one expressed by equation (3.2.11). As equations (3.2.8)-(3.2.10) show the velocity dis- 
tribution becomes highly elongated along the U2 axis, and the shear stress now is inversely 
proportional to the velocity shear, 

«IW~-7^a'. (3.2.13) 

The concept of a viscosity coefficient is not useful in this regime, but if we were to formally 
define u as in equation (3.2.11) we would find that u varies inversely as the square of the 
shear. 

The turnover and reduction of the shear stress when the shear is large is a manifes- 
tation of causality in the theory. From the expression given in (3.2.10) we find that the 
maximum value of |ll]ll2 1 is 

\Wu^\max = 7-7^0-^, (3.2.14) 

2vD 
which shows that the shear stress nuiU2 is always limited to be smaller than the pressure 
(p = na'^/3), regardless of the magnitude of the shear. This automatic stress limiter is 
a very natural and welcome feature of the theory and is an improvement over the simple 
relation (3.2.11). The reason for the reduction of U1U2 at a large shear is not difficult to 
understand. For large values of shear, particles coming from a mean free distance, uit, 
along the xi axis, develop a transverse velocity difference U2 ~ A{uit). Thus for At ^ 1, 
ul ~ t(.|/(Ar)^ ~ (T^/(Ar)^. Therefore the stress, |w]ll2| ~ ufAT ~ a'^/AT, decreases for a 
large shear. 

The asymptotic behavior of the shear stress in the limit of large shear is qualitatively 
very different from that exhibited by the diffusive fiux in particle diffusion. As we have 
shown above, the shear stress actually decreases with increasing shear and goes to zero 
in the limit 2 A ^00. In particle diffusion on the other hand, as the density gradient 
increases the diffusive fiux asymptotes to a constant value which is equal to the product 
of the particle density and the r.m.s. particle speed (Levermore and Pomraning 1981, 
Narayan 1992). 

3.3 Shear Stress from Numerical Simulations 
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The solution written down in equations (3.2.7)-(3.2.10) is valid even for highly non- 
linear steady state shears, where by non-linear we mean that / deviates considerably from 
/o- We have tested this solution using numerical simulations and describe some of our 
results here. 

We set up a shearing cell extending from xi = —0.5 to +0.5 with a large number A^ 
of particles and with boundary conditions imposed so as to enforce the mean velocity to 
follow equation (3.2.1). The particles scatter with a specified mean free time r and are re- 
injected with the appropriate cooling ^ and an isotropic Maxwellian velocity distribution 
in the local rest frame. We let the system evolve until it achieves steady state and then 
compare the results with the theoretical predictions. Figure 4 shows a case with A^ = 
10^, 2A = 1, r = 1, (1 — 0^^ — 1- We display the variation of the moments with time. 
Note that the system achieves a steady state in less than 10 mean free times. 

Figure 5 shows the distribution of the velocity components ui and U2 at the end of 
the run for the three cases that we have simulated. The steady state distribution / is only 
slightly deformed from the Maxwellian post-scattering distribution /o when 2 At = 0.1. 
This is a case of weak shear and it is understandable that / is only slightly perturbed 
from /q. Indeed perturbation theory is valid in this case, and so is the standard weak 
shear formula equation (3.2.11) for the shear stress. In the other two cases that we have 
presented, viz. 2At = 1, 10, we find that / is strongly distorted by the shear. We are 
therefore definitely into the regime of non-linear shear. Nevertheless, we have verified 
that all the second moments calculated analytically (eqs. [3.2.8]-[3.2.10]) are in excellent 
agreement with the numerical simulation. 

We next consider a different situation where the particles have a constant mean free 
path / rather than a constant r. Figure 6 shows UiWi/cr'^ as a function of the shear 
amplitude 2Al/a in this case. The points are results from a simulation of A^ = 10^ particles, 
and the curves were obtained from equation (3.2.10) when we set r = l/a (dashed line) 
or r = 0.55//(T (solid line). Note the interesting result that equation (3.2.10) predicts the 
maximum value of stress accurately for particles with a constant /, although it was derived 
for the constant r case. 

3.4 Steady State Shear with Advection 

The stress-limiter discussed in §3.2 provides a strong indication that our theory satis- 
fies causality. To demonstrate causality in a different and possibly more transparent way 
we consider here a generalization of the previous case where, in addition to the shear in 
V2., we now also allow motion along the xi axis {yi). We refer to this additional motion 
as advection. We do not restrict ^Ji or I}2 in any way except to assume that the gradients 
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of both velocities are in the xi direction and that there are no external forces. Also, we 
retain the assumption of a steady state. The reader may notice a close analogy between 
the flow considered here and that in a steady state accretion disk where xi and X2 would 
refer to the radial [R) and azimuthal {(p) directions respectively. We do not, however, 
discuss rotation in this section (but see §3.6 & 3.7). 

In a steady state, the continuity equation (3.1.5), the X2 component of the momentum 
equation (3.1.6), and the uTu2 component of the second moment equation (3.1.11) give 



1 dn 1 dvi 

ndxi UT^xi 



(3.4.1) 



dv2 _ duiU2 uiU2 dn _ duiU2 uiU2 dvi ('\ao\ 

dxi dxi n dxi dxi vi dxi 

_duTu^ -^dv^ 1 dUi 

Vi^ = -Uf- UiU2-UiU2^ — . (3.4.3) 

OXi OXi T OXi 

Substituting (3.4.2) into (3.4.3) and rearranging terms we find the following expression for 
the steady state shear stress U1U2: 

um = - i ,^,' ) (l-%] nlr^- (3.4.4) 

This relation has several interesting features. 

First notice that if 1)7 = the viscosity coefficient that relates U1U2 to the shear 
91^2/9x1 is just ufr. This agrees exactly with the results of the previous analysis (eqs. 
[3.2.10] and [3.2.8]). When there is advection, the viscosity coefficient is modified by two 
factors. Consider ffist the factor (1 — UT^/tt^), which reduces the viscosity whenever there 
is a steady flow, and is very similar to results previously obtained by Narayan (1992). To 
interpret this result, we need first to ask under what conditions we can have a shear in 
a steady fiow such as we are considering. Since we have ignored external forces, there is 
no outside mechanism that can induce shear in the fiow; the only reason for there to be 
a shear is a downstream boundary condition. For instance, let us suppose that the fiow 
comes in from negative xi with a positive vi and with 1)2 = 0. Let us assume that at some 
large positive xi the fiow meets a sideward moving interface of some sort — a "conveyer 
belt" in the language of Syer and Narayan (1993). The downstream boundary condition 
can be satisfied only if the fiow sets up a shear dWi/dxi. Now, information about this 
downstream requirement has to be transported upstream by the particles themselves. As 
the fiow velocity (as measured in the frame in which d/dt = 0) increases, the net fiux 
of upstream moving particles reduces because only a fraction of the particles has large 
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enough negative velocities to overcome the advection. Therefore, the number of particles 
that participate in the shear stress reduces and this is reflected by a reduction in the 
viscosity coefficient 

Equation (3.4.4) shows that the viscosity coefficient goes to zero when UI = (-u^)^'^. 
In particular, when the shear is weak, we know that -uf ~ cr^/3, and the cutoff occurs at 
Vi = a/v3. The quantity (ufY''^ is the limiting speed of particles in the xi direction in 
our theory. When Ui is equal to this speed, no particles are able to move fast enough to 
beat the advection and as a result no viscous stress is transmitted upstream. This is the 
reason for the cut-off of the viscosity coefficient. In real gases the cut-off may be more 
gradual due to the existence of a high-velocity tail in the particle distribution function. 

We note an apparent paradox in the fact that the viscosity coefficient in equation 
(3.4.4) actually changes sign when vi'^ > uf, which means that the shear stress begins to 
point in the same direction as the velocity shear. It implies that the work done by the shear 
stress, —uiU2dv2/dxi (see equation 3.2.6), is negative, i.e. the shear stress extracts heat 
energy from the flow and converts it into mechanical energy. Since this is not reasonable, 
we interpret the change of sign in (3.4.4) to mean that, when vi'^ > u\ there can be no 
steady state shear in the flow. 

We now come to the second factor (1 + 2Tdvi/dxi)~^ in equation (3.4.4). It shows 
that any divergence in the advection velocity modifies the viscosity coefficient. This is 
fairly straightforward to interpret. Consider first the case when dvl/dxi > 0, which 
corresponds to an expansion of the fiow. In this case, as a particle moves from one region 
of the fiow to another it finds that its velocity becomes closer and closer to the local bulk 
velocity vi. Therefore, the distance that a particle is able to move in "comoving fiuid 
coordinates" before it is scattered is reduced. As a result, the number of particles that 
can reach any point in the fiow from downstream is reduced, and this causes a reduction 
in the shear stress. The case of compression, when dvl/dxi < 0, is similar up to a point. 
In the presence of compression, downstream particles find it easier to move upstream and 
this enhances the shear stress. The interesting point is that the shear stress diverges 
when 2Tdvi/dxi = — 1. What this means is that a fiow with such a large compression is 
unphysical. When 29t'T/9xi < — l/r, the fiow compresses down to a zero volume in less 
than one scattering time. Obviously such a situation can occur only in transient fiows and 
not in a steady state. 

Equation (3.4.4) is a general result for steady fiows which clarifies several physical 
issues. However, because the expression for the stress is given in terms of uf, whose 
magnitude relative to a'^ depends on the strength of the shear, it does not provide a useful 
formula for the viscosity coefficient. To find a more practical formula we restrict ourselves 
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to a steady flow where vi, n and cr^ are constant. In this situation the momentum equation 

gives 

dWu^ = -— ^ (3 4 5) 

dxi dxi ' 

and the four relevant components of the second moment equation (3.1.11) are 

OXi T 6T 

-^ +^7r^ + — + 2W^7r^ = 0, (3.4.7) 



3t 9xi t 9xi 

^ ^^-CT^+W^ + — =0, (3.4.^ 



3r 9x1 



-2'^'^2 , du\U2 , 'U1W2 



WfTT^ + W^T^^ + ^^ = 0. (3.4.9) 

ax\ ax\ T 

From the momentum equation (3.1.6) we find that du\/dxi=Q, and since a^ is in- 
dependent of position, it follows that du^/dxi = —du\/dxi. We now solve equations 
(3.4.5)-(3.4.9) making use of these relations, and find 

which is just equation (3.2.10) modified by the causality factor (1 — 3tJi^/o"^). This provides 
a useful special case of equation (3.4.4) for fiows with uniform advection. 

3.5 Bulk Viscosity 

We now explore the nature of bulk viscosity in our theory. We consider a fiow where 
the mean velocity at each point is oriented along the xi-axis and where all gradients are 
also along xi. Let us define the total time derivative d/dt by 

d d d 

dt = m^''d^,- ('•'•') 

The continuity and momentum equations (ignoring external accelerations) give 

dn dvl ^ /o r r^N 



dvi dp d 

n—— + -; h 



dt dxi dxi 
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n I u 



n 



(3.5.3) 



where in (3.5.3) we have spht the stress tensor into the sum of the pressure p and an 
anisotropic term which wiU be associated with bulk viscosity (cf. eq. [3.1.13]). 

By the symmetry of this problem, there are only three non- vanishing second moments, 
namely uf, u^, u\. Ignoring the acceleration terms, these components evolve according to 

'+ 2-l + -U?---ia2 = 0, (3.5.4) 



dt T 3t 

Adding equations (3.5.4)-(3.5.6) we find 



dt \ dxi tJ St 



M + !i _ l^,^ = 0, (3,5.5) 

dt T 6T 

dUn Wq 1 — f 9 , , 

^ + ^ - -^a^ = 0. (3.5.6) 



—2 dvi f d ^ ^\ 1 2 



^ dxi \dt tJ 2 ^ ' 

This is analogous to equation (3.2.6). The left-hand-side gives the work done per particle 
per unit time. Note that the stress nuf is the sum of the pressure and the bulk viscous 
contribution. The right hand side of (3.5.7) is the rate at which particles acquire kinetic 
energy, some of which is lost through cooling. Equation (3.5.7) thus expresses energy 
conservation. 

Let us first consider the case when the velocity gradient dvl/dxi is small compared 
to 1/r. In this case we expect -uf ~ ti| ~ «§ ~ (t^/3. Equation (3.5.4) gives 



a 



2\ ^T^ Ar..2 



dvi du\ ^ 



ul- — \= -2ul^ _ ^^ _ ^a\ (3.5.^ 

^37 ^dxi dt St ^ 



Substituting u{ ~ a /S into the right hand side and using equation (3.5.7) we find 

n(^-^)--l.on|^, (3.5.9) 

where z/q is the kinematic coefficient of shear viscosity introduced in equation (3.2.12). 
We thus recover the standard result, viz. that in the limit of weak velocity gradients, the 
shear and bulk viscosity coefficients are related by a factor of 4/3. If this relation and 
the analogous equation (3.2.11) are substituted in the momentum equation we obtain the 
Navier-Stokes equation for the evolution of a viscous fiuid. These relations are, however, 
valid only for weak velocity gradients and for slow variations of the fiow parameters. When 
these conditions are violated we need the more complete theory presented in this paper. 
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Let us now proceed to the case when the velocity gradient is not necessarily small. 
For simplicity we assume that the velocity gradient is independent of xi: 

-^ = 2C = constant, (3.5.10) 

ox I 

Let us also assume that the temperature, which is proportional to cr^, is the same ev- 
erywhere and does not change with time. It follows from equation (3.5.7) that uf must 
also be independent of position and time. It is now straightforward to solve equations 
(3.5.5)-(3.5.7), which yields 

Note that ^ is negative for positive C. This is because particles cool when there is an 
expansion of the flow and we need a source of heating (rather than cooling) to maintain 
the moments at a constant value as we have assumed. 

For a small velocity gradient, \Ct\ <^ 1, equation (3.5.12) recovers the usual result 

^- — ^-- — -2(7, (3.5.13) 

^333' ^ ^ 

which is identical to equation (3.5.9). When Ct is large and positive we flnd ^ -^ —1/2 
and Ui/a'^ -^ 0. For an extremely rapid expansion the particles are very cold in the xi 
direction and their velocity distribution is restricted nearly to the U2-Us plane. This is 
reasonable. When Ct is large and negative, equations (3.5.11), (3.5.12) appear to reveal 
a problem. If Ct < —1/4, then we flnd uf > a'^ which is physically inconsistent. There 
is, however, a simple explanation for this result. Precisely when Ct = —1/4 we notice 
that ^ = 1. Going back to the original equation (3.1.4) where ^ was deflned, we see that 
^ = 1 is a physical upper limit to the cooling, representing the case when all the energy 
in a particle is removed entirely in a collision. However, Ct = —1/4 corresponds to an 
inflnite compression of the gas in one collision time, which generates an inflnite amount of 
heat. Clearly our cooling term is unable to remove this heat, and thus the assumption of 
a steady state breaks down. The present limit on Ct is exactly equivalent to the limit on 
9^^/9x1 discussed in §3.4, and these have the same underlying physics. 

3.6 Shear Viscosity in DiflFerential Rotation 

We consider flnally a shearing rotating flow as an example of an application with 
velocity-dependent accelerations. This is also of particular importance in astrophysics 
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because of applications in accretion disks. Kato and Inagaki (1993) discuss rotating flows 
in greater detail than we do, but they have limited themselves to a small shear, whereas the 
following discussion is applicable even for a large shear. (See also Goldreich & Tremaine 
1978, who have derived an expression for the shear stress in a rotating disk of particles in 
Saturn's ring using a different approach.) 

We consider a rotating flow in a local Cartesian approximation, called the shearing 
sheet approximation (Goldreich & Lynden-Bell 1965, Narayan, Goldreich & Goodman 
1987). Let the flow take place in a central potential which produces an inward radial 
acceleration g{R). Consider a reference radius Rq and let O = [g{Ro)/Ro]^^'^ be the 
equilibrium angular velocity such that centrifugal acceleration just cancels g aX R = Rq. 
We work in a frame which rotates with angular velocity O and is centered on Rq. We 
take our local Cartesian grid to have xi = R — Rq, X2 along the azimuthal direction, and 
xs parallel to O. We have arranged so that a particle at rest at xi = X2 = ^3 = is in 
equilibrium. However, particles at non-zero xi experience an additional acceleration of the 
form g'xi, arising from a combination of the central potential and centrifugal force; we 
have g' = dg/dR + O^. In addition, moving particles in the rotating frame experience a 
Coriolis force. Therefore, the three components of the acceleration are given by 

ai = g'xi + 20i;2, «2 = —20^1, as = 0. (3.6.1) 

For a truly spherical central potential, a^ 7^ 0, but we have simplifled matters by assuming 
a cylindrically central potential (as appropriate for a thin disk). 

Let us consider a steady state shear flow in this rotating frame, and let us assume 
that n is constant and that (cf. §3.2) 

IJY = 0, v^ = 2Axi, W = 0. (3.6.2) 

The four non- vanishing second moments, uf, u^, u^, and U1U2, satisfy the following equa- 
tions (compare to equations 3.2.2 — 3.2.5): 

3wf-(l-O^^-12OrIITII^ = 0, (3.6.3) 

3^-(l-0^^ + 12SrWw2 = 0, (3.6.4) 

3^ -(1-0^^ = 0, (3.6.5) 

2Srwf-20rw| + 11111^ = 0, (3.6.6) 
where we have deflned the vorticity 2B by 

2B = 2A + 2n. (3.6.7) 
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As before, summing equations (3.6.3), (3.6.4) and (3.6.5) we obtain the usual energy con- 
servation relation 

£ 1 

-Wu^ ■ 2A = - ■ -a^ . (3.6.8) 

r 2 

Solving the equations for the cooling constant ^ and the various second moments we find 

^ " 3 + 8A2r2 + 12/^2^2' ^^-^-^^ 

- l + 80V + 2/.V^ 2 

"^ = 3 + 8^2.2 + 12,2,2 ^ ' (3-6.10) 

- 1 + 8i?V^ + 2^2,2 

^^ = 3 + 8A2r2 + 12^2,2 ^ ' (3-6.11) 

2Ar 2 ^ X 

^^"^ ^ -3 + 8A2r2 + 12^2,2^ ' (3-6-13) 

where k^ is the epicyclic frequency, 

k'^ = AQB. (3.6.14) 

These are generalizations of the results of §3.2 when there is rotation. We notice 
that the results now depend on two dimensionless parameters. At and kt, rather than 
one, because the flow now has two frequencies associated with it, the shear 2A and the 
epicyclic frequency k. 

We may consider various limits of equation (3.6.13). If Ar, kt ^ 1, which corresponds 
to very rapid collisions, we recover the usual result 

uTu^= — ■2A = -uq—^, (3.6.15) 

3 oxi 

which is the standard formula for shear viscosity. In this limit, rotation plays no role. 
Similarly, if kt <^ 1 and if At ^ 1, i.e. if we have strong shear in the presence of weak 
rotation, we again recover equation (3.2.10). On the other hand, if kt ^ 1, then even if 
the shear is weak (i.e. At <^ 1), the shear stress is suppressed relative to the non-rotating 
case (see also Kato and Inagaki 1993), 

^^^ ^ - 1^] 2 2 ■ V ■ 2^- ^3.6.16) 

1 + 4K^r^ 3 

The reason for this is straightforward. When k ^ 1/r a particle undergoes many epicycles 
within a mean free time and travels a distance equal only to the radius of an epicycle in 
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the Xi direction. Therefore, the effective mean free path is much shorter than that in the 
non-rotating case and this suppresses the viscosity. FinaUy, if At, kt are both ^ 1, then 
the viscous stress is reduced by an even larger factor. 

In the derivation of the basic equations we had set the third moment of the velocity 
to zero. Physically this condition is satisfied when the velocity distribution of particles has 
a reflection symmetry, or in particular when the distribution is constant on ellipses. This 
turns out to be a very good approximation for shear flows in rotating systems, where the 
epicyclic motion of particles results in a quite symmetric velocity distributions (see figure 
8). 

Obviously, the shear stress nuiU2 in the presence of rotation is limited by the two 
parameters discussed above, instead of just one as in the non-rotating case (eq. [3.2.14]). 
The maximum stress is thus smaller than in non-rotating flows. We discuss the stress 
limiter further in §3.7, where we specialize some of these results to the case of a Keplerian 
disk. 

To demonstrate causality more explicitly, we restore the temporal and spatial deriva- 
tives in the V2 component of the momentum equation and the equation for U1U2: 

dv^ d 



dt dxi 



{Wu^), (3.6.17) 



— ^- ^^1-Z \- [uf - u^)2il = — U1U2. (3.6.18) 

We continue to assume that n is constant and that t^T = 0. This assumption is technically 
inconsistent since variations in V2 will cause fluctuations in vi (through the Coriolis accel- 
eration) which will lead to compressive motions and changes in n. However, compression 
in n will cause sound waves which will complicate the analysis. We fllter out the sound 
waves by ignoring variations in vi and concentrating only on V2 and uiU2- Let us further 
assume that the velocity shear ^I'i'/^xi is small. In this limit, we have 2At <^ 1, S ^ O, 
K ~ 20. Equations (3.6.10) and (3.6.11) then indicate that u^ ^ U2 ^ a'^/3. Substituting 
these relations into (3.6.18), we flnd 

We now consider small perturbations in V2 and uTu2' 

V2^V2+SV2, U1U2 — > U1U2 + duiU2- (3.6.20) 

Equations (3.6.17) and (3.6.19) then combine to give 



dt 3 dx'^ a"^ dt'^ 
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(3.6.19) 



This is analogous to equation (2.3.11) for particle diffusion and has the familiar wave 
operator which enforces causality. Thus, we have proved that even in the presence of 
strong rotation, viscous signals travel at a finite speed and satisfy causality. 

3.7 Limits on a— Viscosity 

Shakura and Sunyaev (1973) used simple scaling arguments to write the kinematic 
viscosity coefficient in a turbulent thin accretion disk in the form 

v = a--^, a<l, (3.7.1) 

\Ik 

where Q.k is the Keplerian angular velocity. Their argument was that viscosity depends 
on two quantities: (i) the speed a of the turbulent eddies, which is likely to be limited by 
the sound speed c^, and (ii) the mean free path / which is likely to be no larger than the 
vertical scale height of the disk, the latter quantity being ~ Cg/Qx- Since u ~ cr/, they 
thus obtained equation (3.7.1) with a limited to be < 1. 

In §3.6 we have developed a fairly complete treatment of viscous interactions between 
particles in a differentially rotating system and we can therefore quantify equation (3.7.1) 
more carefully. Our theory has two parameters. One of these is the r.m.s. particle velocity 
cr. In the context of disk viscosity, our "particles" are turbulent eddies or blobs, and in a 
purely hydrodynamic situation we expect the blob velocity a to be smaller than the local 
sound speed, otherwise it would be dissipated by shocks. Let us therefore write 

a = ttiCs, tti < 1. (3.7.2) 

Our second parameter is the mean free time r between collisions. Since we cannot a priori 
assign any limits to r, let us simply express it in terms of the rotation frequency O, 

r = f , (3,7.3) 

where we do not restrict a2 to any particular range. We showed in §3.6 that the viscosity 
depends on two frequencies, A and k, which are given by 

2A Rdn dlnn , „ ^, 

— = = f3 7 4l 
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Substituting these relations into equation (3.6.13) we find the following result for the 
coefficient of viscosity 

^^^ ^ = 0-^ ^3 7 61 

3 + 2al{d\nn/dlnR)'^ + 24al{dlnnR'^/d\nR)n ~ O' ^ ' ' ' 

which gives an expression for the Shakura-Sunyaev a parameter in terms of cti and a2- In 
the particular case of a Keplerian disk, where fl{R) oc i?~^/^, we have 2A/fl = —3/2 and 
K? /Q? = 1, and (3.7.6) becomes 

UK = TTl- TT = CtKTT- 3.7.7 

The most interesting feature of these relations is the fact that u and uk each has an 
absolute maximum value regardless of the value of a2. Equation (3.7.7) for instance reaches 
its maximum value when Or = a2 = a/2/11 at which point we have 

2 

{aK)ma. = -^ = 0.071a?. (3.7.8) 

This is an extremely small limit for ax considering that cti is itself limited to lie below 
unity. Note that if a is defined alternatively through the stress tensor relation uiU2 = 
—acl, then the above bound should be multiplied by 3/2. 

To demonstrate how small aK is in a Keplerian disk we show in figure 9 the variation 
of ax as a function of the parameter a2 = Or. We see that, if a2 lies more than a factor 



of few outside the optimum value of a/2/11, then aK falls even further. Also, cti itself 
is presumably somewhat less than unity since turbulent blobs are likely to have a range 
of speeds with the fastest blobs limited to speeds less than Cg. Thus, we conclude that 
turbulent viscosity in accretion disks, when expressed in the Shakura-Sunyaev form (3.7.1), 
is likely to have values of uk which are somewhat smaller than 0.07. Lower values of uk 
imply disks with larger surface densities and spectra that are more nearly thermal. 

In regions away from the Keplerian region of a disk, for example in the inner boundary 
layer, we need to use the full expression given in equation (3.7.6). We see that in regions of 
large shear the effective a will be much smaller than in the Keplerian case. Figure 9 shows 
a case where d\nVt/d\nR = 5, a modest value by the standards of boundary layers (e.g. 
Narayan and Popham 1993, Popham et al. 1993). Note how small the viscosity coefficient 
becomes in this case. When we add to this the fact that the boundary layer region has 
radial flows which also reduce u through the causality factor discussed in §3.4, we see that 
accretion disk boundary layers will in general have very weak viscosity. The effect that 
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this will have on the structure of the boundary layer is yet to be investigated. The role of 
third moments, which may not be negligible in boundary layers, is also unclear. 

The formulae we have developed are based on the assumption of a constant mean free 
time, and one may wonder whether the results would change significantly with other scat- 
tering laws. However, to leading order, different scattering laws should all be describable 
through an effective mean free time, and since we did not restrict Or in any way our limit 
on ot would appear to be robust. The point is that, ultimately, the limit on a comes from 
the stress-limiter, which inhibits the viscous stress from exceeding the pressure, coupled 
with geometrical factors which further limit the shear stress because of the epicyclic mo- 
tion. These are sufficiently general principles that we cannot envisage violating the limits 
shown in figure 9 by any significant factor. Obviously, if the shear stress is produced by 
non-hydrodynamic effects, e.g. magnetic fields, then our limit could be exceeded. 
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4. SUMMARY AND DISCUSSION 

In this work we have attempted to incorporate causahty into the theory of transport 
phenomena. We have developed equations for particle diffusion in one and three dimensions 
which are an improvement over the standard diffusion equation whenever gradients are 
large and causality limits the flux. We have also generalized the equations of viscous 
hydrodynamics to include situations where velocity gradients are large or time variations 
are rapid. In such cases the standard Navier-Stokes equation breaks down. 

Our approach is based on taking moments of the Boltzmann equation and using an 
approximate form for the collision term (equation 2.1.3). We assume that particles collide 
with a mean free time r which is independent of velocity. Furthermore, we assume that 
each scattering completely isotropizes the velocity distribution of the particles. These 
assumptions permit us to write simple expressions for the velocity moments of the post- 
scattering distribution function /o in terms of the moments of the pre-scattering function 
/. The theory is constructed to handle situations where / differs appreciably from /q. 
This is an important departure from the usual perturbative method of analyzing transport 
phenomena where the quantity (/ — /o) in equation (2.1.3) is assumed to be small compared 
to the equilibrium distribution /q. The non-perturbative nature of our theory allows us 
to write down moments of the Boltzmann equation which are valid even in situations that 
involve large departures from equilibrium. 

The first example we have presented is the case of light particles diffusing in one 
dimension in a fixed background (§2.1, §2.2). We take the zeroth and first moments of 
the Boltzmann equation and close the equations by setting the second moment of the 
velocity, v^ ^ to a constant. This leads to a modified diffusion equation (2.1.11) which 
differs from the usual diffusion equation in having a wave operator instead of the standard 
Laplacian operator. The presence of the second time derivative in the wave operator 
enforces causality (Israel and Stewart 1980, Schweizer 1984). This term arises from a time 
derivative term in the first moment equation, which is neglected in deriving the standard 
diffusion eqution. We therefore deduce that, if we wish to preserve causality, it is important 
to retain consistently all time derivatives of higher moments. We follow this prescription 
in the various problems we have analyzed and we find that we obtain equations that are 
causal in all cases. 

The Green's function of the one dimensional diffusion problem reveals several inter- 
esting features. First, it vanishes outside the particle "velocity cone," \x\ > at, where 
X is distance from the point of injection of the particles and a is the r.m.s. velocity of 
the particles. This explicitly demonstrates the causal nature of the equations. Secondly, 
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the Green's function has J-functions exactly at the two edges of the velocity cone. This 
shows that the theory effectively replaces the full velocity space distribution f{v) by two 
5-functions at v = ±a. This simplification means that the theory cannot treat phenomena 
like phase mixing which arise from a continuous distribution of velocities. If such effects 
are important in any application then the theory we have presented is inadequate. One 
should either integrate the single velocity Green's function over the particle distribution 
function or use higher order moments (see equations [2.1.19] & [2.1.20]). The theory also 
ignores the tail of high velocity particles which are normally found in a real distribution 
function f{v). Real gases do not diffuse with an exact cutoff in space. While the bulk of 
the signal is limited to a speed equal to the r.m.s. speed a of the particles, there is always 
an exponentially small number of particles at large velocities propagating faster than a. 
By eliminating the high velocity tail altogether our theory ignores this weak precursor 
signal and therefore has an exact cutoff of signals at a maximum propagation speed a. 

Leaving aside the issues of phase mixing and the exponential tail, we show that the 
theory does well in describing the behavior of the bulk of the particles. In fact, when 
all particles have a constant mean free time r the theory gives the mean square distance 
{x^) traversed by particles exactly, regardless of their velocity distribution function. The 
agreement is perfect both at early times, when the particles are streaming freely, and at 
late times, when they have undergone many scatterings. Furthermore, even for a more 
general scattering law we can usually find an equivalent r such that our theory provides a 
satisfactory description of the true evolution of (x^) (see figure 2). 

It is a straightforward matter to generalize the theory to diffusion in three dimen- 
sions. The only complication originates from the fact that the minimal extension of the 
one-dimensional closure relation requires a somewhat more extreme assumption in three 
dimensions, namely that the velocity second moment tensor is not only constant but also 
isotropic. This assumption leads to a causal but somewhat unphysical Green's function 
with negative densities. Despite this unsatisfactory feature, the theory provides an excel- 
lent description of the evolution of the mean square distance (r^) travelled by the particles. 

An interesting special case of the three dimensional problem is when the density gra- 
dient is limited to a single direction, say the xi axis. In this case, we find that our three- 
dimensional diffusion equation reduces to a form that is identical to the one dimensional 
equation except that the speed of propagation of signals is no longer the r.m.s. particle 
speed a but rather a/v3. In this particular case, it is clear that we effectively replace the 
full distribution function /(t'l, ^2, ^3) by one that is restricted to two planes in velocity 
space, i.e. vi = ±(t/v3- We have already seen the limitations associated with such an 
approximation, namely the inability to model phase mixing and signal propagation by the 
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exponential tail of fast particles. But apart from this, we expect the one-dimensional pro- 
jection of the three-dimensional equation to behave physically and to predict the behavior 
of the majority of particles in a satisfactory way. In particular, the Green's function for 
this problem does not have unphysical negative densities. 

In §3 we have used our method to derive a general causal equation (3.1.11) for the 
evolution of the stress tensor. Since the stress tensor is of second order in the velocities, 
we take the zeroth, first and second moments of the Boltzmann equation, assuming as 
in the diffusion problem a constant mean free time r. We also set the third moments of 
the velocity to zero, which allows us to close the moment sequence and to derive a self- 
contained set of equations. The main additional feature in our model of the scattering is 
the introduction of a "cooling" factor ^ (see equation 3.1.4) whereby scattering events are 
allowed to be inelastic. This modification allows us to model cooling processes where some 
of the kinetic energy of colliding particles is converted into another form of energy, e.g. 
radiation, and removed from the system. By including cooling, we are able to investigate 
steady state problems even when there is heating due to stresses. 

Before proceeding to describe our results on the stress tensor, we discuss the effect 
of the various approximations we have made. The physical systems we have investigated 
involve velocity gradients along a single direction, a case where we can expect the theory 
to be well-behaved since we know that the particle diffusion problem admits a physical 
solution. The theory assumes a constant mean free time for all particles and therefore 
involves a simplification of the real collision integral. However, as we have shown in 
detail in the diffusion problem, and with some limited numerical experiments in the stress 
problem as well (figure 6), general scattering laws can usually be reduced to an effective 
constant r with very little qualitative differences. Therefore, the assumption of a constant 
r is unlikely to be a serious limitation. 

The most serious simplification of our theory is the neglect of third moments. This 
approximation is reasonable when the velocity distribution has a refiection symmetry, or 
when the mean velocity profile is nearly linear in the vicinity of the region of interest. As 
we show in the Appendix, the third moment terms will be small provided the dimensionless 
second velocity derivative, ar'^d'^v/dx^, is small compared to either unity or the square of 
the dimensionless shear (equation AlO). There is another similar condition on the third 
derivative of the velocity (equation All). These conditions may not be satisfied in highly 
transient fiows or in the vicinity of shocks, but should be valid in smooth steady state fiows 
and especially in rotating systems where the distortion of the particle velocity distribution 
tends to be limited (see figure 8) . It is important to emphasize that the neglect of the third 
moments does not require the velocity gradient to be small; the theory we have developed 
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is valid for arbitrarily large gradients as long as the gradient itself is spatially constant, 
i.e. the second derivative of the velocity is sufficiently small. 

Subject to the above caveats, equations (3.1.6) & (3.1.11) may serve as potential 
candidates to replace the Navier-Stokes equation whenever spatial or temporal derivatives 
are large. We have not tested these equations in their full generality. Rather, we have 
concentrated on special cases involving steady state shear flows to demonstrate several 
physical effects involving causality in viscous transport processes. 

In §3.2-3.4, we derive some general results for the shear stress in a steady shear flow. 
According to the standard theory, a velocity shear dv^jdxx = 2A gives rise to a shear 
stress nu{2A), where n is the density and u is the kinematic coefficient of viscosity. Our 
causal theory reproduces this result in limit that A is small, with z/ = a^r/d, but more 
generally reveals three distinct new effects described below: 

1. In §3.2 we consider steady shear of arbitrary magnitude 2A without advection and 
show that the shear stress is given by nu{2A){l + 2[2Ar]^/3)~^ (equation 3.2.10). 
This expression gives a significantly smaller stress than the standard result mentioned 
above when the dimensionless shear 2At is large. Whereas the standard relation for 
the shear stress diverges in the limit of large shear, the modified expression reaches a 
maximum value, which is less than the pressure, at a finite value of the dimensionless 
shear 2At and reduces for larger values of the shear. This behavior is a consequence 
of causality. The particles at any given point typically have originated a mean free 
path away. Therefore, for a large shear, they have a large transverse velocity spread 
w| ~ 0"^ ^ ul- This translates to a stress limit |wTw2| < {uf v3^^l'^ <^ a^. We have 
tested our result for the shear stress through numerical particle simulations spanning 
a range of values of the dimensionless shear 2At (§3.3). The agreement is essentially 
perfect (to within numerical precision) when all particles have the same mean free 
time. Even when we instead consider the mean free path to be constant, our theory 
agrees well with the numerical simulations (figure 6) . In particular, the theory predicts 
the maximum value of the shear stress accurately. Note that in the case of the shear 
stress, not only does causality set a limit to the maximum stress, it actually causes a 
reduction of the stress in the limit of large shear. This asymptotic behavior is very 
different from the behavior of scalar particle diffusion (e.g. Levermore and Pomraning 
1981, Narayan 1992), where the fiux asymptotically approaches to a constant value in 
the limit of a large density gradient. 

2. When a constant advection velocity, Ui, is added to the steady shear flow, we flnd 
that the shear stress is reduced further by a second suppression factor (1 — 3vi'^ /uf) 
(cf. eq. [3.4.4]). When there is advection, the shear stress has to be communicated 
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by particles that move upstream faster than the advection speed Ui. The flux of such 
upstream moving particles reduces as the advection speed increases and cuts off when 
vi exceeds the r.m.s. velocity (ttf)^/^ of particles in the direction of the flow. A 
similar qualitative result was obtained by Narayan (1992) who showed that causality 
leads to a reduction and ultimate cutoff of viscosity when there is advection. Narayan 
considered somewhat general velocity distribution functions whereas our discussion 
corresponds to a very simple velocity distribution consisting of two 5— functions at 
ui = ±(tt^)^/^. On the other hand, Narayan assumed that the shear is small, whereas 
we allow arbitrarily large velocity gradients. One other difference is that Narayan's 
expression for the viscosity coefficient left undetermined the critical advection velocity 
at which the viscosity will vanish. The present theory shows that, at this level of 
approximation, the cutoff occurs exactly at the sound speed a j \pi. We caution, 
however, that in real gases the advection factor may not involve a sudden truncation 
of the viscosity coefficient but rather an exponential cut off as IJx exceeds (wf )^' ^. This 
is because of the high velocity tail in the particle distribution function which may be 
able to transmit a weak signal upstream even at large advection speeds. 
3. We have also investigated the effect of a non-uniform advection and find that the shear 
stress is modifled by a third factor, (1 ^2Tdv\ldx\)~^ (cf. eq. [3.4.4]). This shows 
that the shear stress is suppressed when the flow expands, i.e. when dvxjdxx > 0, 
but is enhanced in the presence of compression. Curiously, the stress diverges when 
2Tdvi/dxi = —1, but this is not a serious problem since it refers to an unphysical 
limit where the gas is compressed to an inflnite density in a collision time. The same 
limiting compression appears also in the context of our analysis of bulk viscosity in 
§3.5. Kato and Inagaki (1993) have independently noted that expansion or contraction 
of the advecting flow can modify the shear stress. 

In §3.6, we extend the analysis to rotating flows and flnd that the above results for 
the shear stress are further modifled. We flnd that the shear stress in a rotating flow in the 
absence of advection is given by ?iz/(2A)(l + 2[2Ar]^/3 + 4[Kr]^)~^, where O is the angular 
velocity and k'^ = 40(A + 0) is the square of the epicyclic frequency. Thus, in the presence 
of rotation, the shear stress is reduced compared to the classical diffusion formula nu{2A) 
when either the shear or the epicyclic frequency is large compared to the collision frequency 
1/r. The suppression effect due to k has been noted in previous studies by Goldreich and 
Tremaine (1978) and Kato and Inagaki (1993). Kato and Inagaki's approach is very similar 
to ours, whereas Goldreich and Tremaine treat two particle inelastic collisions in much 
greater detail, but they close the moment equations by assuming that the distribution 
function in velocity space has a triaxial Gaussian shape. When kt ^ 1, a typical particle 
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undergoes many epicycles within a collision time and consequently moves radially only 
across the limited scale of an epicycle. The effective mean free path is therefore smaller 
than in the non-rotating case by a factor ~ kt, and the stress is reduced by the square of 
this factor. For the same reason the distribution of particle velocities in phase space does 
not stretch indefinitely as in the non-rotating case (figures 5 and 7), but rather reaches a 
limiting oval shape as O increases beyond 1/r (figure 8). The new feature that comes out 
of our theory is that the viscosity coefficient is suppressed not just by the dimensionless 
epicyclic frequency kt but by a combination of the dimensionless shear 2 At and kt. This 
result is of potential importance in astrophysics. For instance, in Keplerian disks, the two 
frequencies 2A and k are of comparable magnitudes. 

Finally, in §3.7, we have applied our results to thin accretion disks, where we find 
that the causal limit on the shear stress leads to a constraint on the value of the disk 
viscosity coefficient u. This limit can be translated to an upper bound on the dimension- 
less a-parameter defined by z/ = ac^/O, where Cg is the sound speed in the disk. For 
hydrodynamic turbulence in a Keplerian disk we find that a < af/ylOS < 0.07 for any 
value of Or (see figure 9), where cti < 1 is the ratio between the r.m.s. turbulent blob 
speed and the sound speed in the disk. This bound can be exceeded only if the turbulence 
is non- hydrodynamic (e.g. magnetic) and cti > 1. Note that our upper bound should 
be scaled up by a factor of 3/2 if a is instead defined through the stress tensor relation 
U\U2 = ^otc^s- "^^^ upper limit on a is much stronger in boundary layers where the velocity 
shear is larger than Keplerian. 

To conclude, we note three issues related to viscous flows that we have not addressed 
in this paper, namely (i) the effect of advection in rotating flows (where the physics appears 
to be a little different than in the non- rotating case because of the Coriolis acceleration), 
(ii) the structure of viscous shocks, and (iii) time-dependent effects. However, the basic 
equations of our theory, especially the stress-evolution equation (2.1.11), appear to be 
capable of dealing with these issues. Work along these directions may provide further 
insight into the role of causality in astrophysical accretion disks. 
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Appendix 

In the following we rederive tiie results given in §3.2 for the components of the stress 
tensor in the presence of an arbitrary velocity shear, dWijdxx = 2A. We employ here a 
different, and possibly more transparent, approach in which we follow the trajectories of 
individual particles. We then extend this approach to investigate under which conditions 
our assumption of vanishing third moments is valid. This helps us identify the limit of 
validity of the theory developed in the paper. 

Consider the particles at a particular plane, xi = 0, at a given time and label the 
particles by the coordinate xi at which they last suffered a scattering. Let the relative 
velocity component of a particle immediately after the last scattering be given by ('Ui)o, 
{u2)q, {us)o- The relative velocity components at position xi = are then 



2 aa^ 



ui = {ui)o, U2 = {u2)o + 2Axi + --^-^xj---, Us = {u3)o, (Al) 



where we have Taylor-expanded the mean velocity V2 as a function of xi. By the assump- 
tion that the scattering is isotropic, we have 



(uDo = {uDo = (uDo. (A2) 

Furthermore, since ui and us of a particle do not vary in between scatterings, we have 



ui = K)o, ui = iui)o. (A3) 

Consider a particle at xi = with a particular value of ui. Given our assumptions of 
constant density, constant second moments, and constant mean free time, the probability 
that this particle had its last scattering between xi and xi + dxi is given by 

p{xi)dxi = - — — exp{xi/uiT)dxij xiUi < 0. (A4) 

\ui\t 

We then see that, averaged over all particles at xi = 0, odd moments of xi such as (a;i) 
and (a;f ) vanish, while even moments become 

{xD = 2^y, {xj) = 24^r^ ■ ■ ■ (AS) 

Let us consider now the case discussed in §3.2, where the velocity profile is exactly 
linear with a slope 2A and where the second derivative of the velocity is zero. From 
equation (Al), we find the second moment of U2 to be given by 



ul = {ul)o + 4A''{x'i) = (1 + 8A^T^)ul (A6) 
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Using now the relation, a^ = w^ + w| + tig, we find that 



a 



m 



ui 



Ui 



1 + 8AV^ 
3 + 8A2r2 ' ^^ 3 + 8A2r2 

Similarly, we find the velocity moment U1U2 to be given by 



a 



U1U2 = {uiU2)o + {2Axi{ui)q) = -2Atu\ 



2At 



:Cr 



(AG) 



(A7) 



3 + 8A2r2 

These relations are identical to the results derived in §3.2. 

Let us proceed next to a case where the second velocity derivative d'^V2/dx\ and 
higher derivatives are not zero but small and let us estimate the contribution that they 
make to the shear-stress. This is best accomplished by considering the U1-U2 component 
of equation (3.1.11) which is given below: 



^ .^j du^U2 

U1U2 = —IrAui — T—-^ — , 

ox I 



(AS) 



We can estimate U1U2 using equations (Al) and (A4). Substituting this into equation (A8) 
we find that the modification to the shear stress is small if the following condition on the 
derivatives of V2 is satisfied: 



d 



dxi 



' 2— ^^^2' 
r u^ 



^ dxl 



< 2Aul, 



(A9) 



Using the Gaussian closure relation, uj = 3('uf) , this means that we require 



and 



ar 



a\^ 



d'^V2 



dx\ 



2^2^ 



< 



(3 + 9>A't 



d^v^ 



dx\ 



2^2^ 



< 



2Ar(3 + 8A2r 



(AlO) 



(All) 



We thus see that, when the dimensionless shear 2 At <^ 1, the condition for the third 
moments to have a negligible effect is that (i) the dimensionless second velocity derivative 
in the left-hand side of equation (AlO) must be small compared to unity, and (ii) the 
dimensionless third derivative in (All) must be small compared to 2At. When 2At ^ 1, 
both of these criteria are weakened by a factor of A^r'^. 

For a rotating system the value of xi is limited by the radial extent of an epicycle 
~ -ui/k. In the case of a thin disk with a scale- height h <^ r and a Keplerian velocity 
profile, the conditions (AlO) and (All) will be automatically satisfied since the left hand 
side is always of order {h/r) ^ 1. Thus the limit on a which was derived in §3.7 under 
the assumption that the third velocity moments vanish is self-consistent. 
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FIGURE CAPTIONS 

Fig. 1: Comparison between tiie Green's function for the one- dimensional diffusion of a 
Maxwellian distribution function (dot-dashed line), the Green's function of the standard 
diffusion equation (dotted), and the Green's function (2.1.13) of the modified diffusion 
equation (2.1.11) (sohd). The curves are normalized to a unit area and are shown at two 
different times: t = O.lr (upper panel) and t = lOr (lower panel), where r = const is 
the collision mean free time. In the lower panel, the two 5— functions which propagate at 
V = ±a have not been plotted since they are suppressed by e"*/^"^ = e~^. 

Fig. 2: The evolution of the r.m.s. particle distance (x^)^/^ for one-dimensional diffusion 
away from an instantaneous point source. Length and time are normalized by the mean 
free path / and the mean free time //a, respectively. The different curves compare the exact 
result for a one-dimensional Maxwellian distribution of particles with a constant mean free 
path (solid line) to our one-dimensional causal model that assumes a constant mean free 
time r = (2/7r)^/^//(j (short dashed) and the standard diffusion approximation result (long 
dashed). The causal model fits the true behavior very well. Note that our model describes 
the evolution of (x^) precisely for any system of particles that has a constant mean free 
time. 

Fig. 3: The same as in figure 2 but for the r.m.s. particle distance (r^) in three- 
dimensional diffusion. 

Fig. 4: Evolution of the second velocity moments with time. The results were obtained 
from a numerical simulation of a linear shear fiow with 10^ particles, corresponding to 
2A = 1, (1 — 0^^ = 1? ciiid r = 1. Note that the moments achieve their steady-state values 
uf = u"^ = 1/3, w| = 1 and uiU2 = —1/3 in less than 10 mean free times. 

Fig. 5: Phase-space distribution of particles in the U1-U2 plane for steady-state fiow with 
a linear shear. Results are shown for three different shear amplitudes 2At, namely: (a) 
0.1, (b) 1.0, and (c) 10, where r is the mean free time, taken to be independent of particle 
velocity. Note the strong distortion of the velocity distribution at large shear. 

Fig. 6: The stress in a steady-state shear fiow of particles with a constant mean free 
path /. The points show uiv^/cr'^ as a function of the shear amplitude 2Al/a, according to 
numerical simulations with 10"^ particles each. The curves show the prediction of equation 
(3.2.10) when we set r = l/a (dashed line) or r = 0.551 /a (solid line). The maximum 
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stress is predicted accurately by our model even though it assumes a constant mean free 
time T instead of a constant /. 

Fig. 7: Phase-space distribution of particles in the U1-U2 plane for steady-state flow of 
particles with a linear shear 2 A and a constant mean free path /. The different panels 
correspond to different shear amplitudes 2Al/a, namely: (a) 0.1, (b) 1.0, and (c) 10. 

Fig. 8: Phase space distribution of particles in the U1-U2 plane for a steady shear in a 
rotating Keplerian disk. Results are shown for four values of fir, namely: (a) 0.1, (b) 
0.4264, (c) 2.0, and (d) 10. Case (b) yields the maximum value of uTu2- Note the limited 
distortion of the velocity distribution for large shear, in contrast to the cases shown in 
figures (5) and (7). This is because the asymptotic form of the velocity distribution in 
a rotating flow under a large shear is moderated by the epicyclic motion (see equations 
[3.6.10]-[3.6.12]). 

Fig. 9: The value of the a-viscosity parameter for hydrodynamic turbulence in thin disks 
(cf. eq. [3.7.6]) as a function of 02 = ^t, where O is the disk angular velocity and r is the 
mean time between collisions of turbulent blobs. The value of a is divided by af , where 
cti = a/cs < 1 is the ratio between the r.m.s. blob velocity and the sound speed in the 
disk. Results are presented for two rotation proflles: a Keplerian disk (solid line) and a 
"boundary layer" profile with 2 A — Sfi (dashed). 



44 



